clear
use gender_cty.dta

gen cty_fips=cty


merge 1:m cty_fips using causal_normalasylum.dta
rename _merge mergefamstructure

#delimit ;
local vars w2_pos_30_q5_m_tp w2_pos_30_q5_f_tp w2_pos_30_q4_m_tp 
w2_pos_30_q4_f_tp w2_pos_30_q3_m_tp w2_pos_30_q3_f_tp w2_pos_30_q2_m_tp 
w2_pos_30_q2_f_tp w2_pos_30_q1_m_tp w2_pos_30_q1_f_tp 
w2_pos_30_q5_m_sp w2_pos_30_q5_f_sp w2_pos_30_q4_m_sp 
w2_pos_30_q4_f_sp w2_pos_30_q3_m_sp w2_pos_30_q3_f_sp w2_pos_30_q2_m_sp 
w2_pos_30_q2_f_sp w2_pos_30_q1_m_sp w2_pos_30_q1_f_sp;

#delimit cr
forvalues i = 1(1)5{
local gen m f
foreach z of local gen{
local par tp sp
foreach x of local par{
gen est`i'_`z'_`x' = .
gen lb`i'_`z'_`x' = .
gen ub`i'_`z'_`x' = .
gen meanasylum`i'_`z'_`x' = . 

#delimit ;
local vars w2_pos_30_q`i'_`z'_`x';

foreach w of local vars{;
di "`w'";

reghdfe `w' hasnormalschool if hasnormalorasylum==1, absorb(statechetty) cluster(statechetty);
#delimit cr
 mat A = e(b)
 mat B = e(V)
gen df = e(df_r)
gen se = sqrt(B[1,1])
replace est`i'_`z'_`x' = A[1,1]
replace lb`i'_`z'_`x' = est`i'_`z'_`x' - (invt(df,.975)*se)
replace ub`i'_`z'_`x' = est`i'_`z'_`x' + (invt(df,.975)*se)
su est`i'_`z'_`x' lb`i'_`z'_`x' ub`i'_`z'_`x'
matrix drop A B 
drop df se
su `w' if hasnormalorasylum==1 & hasnormalschool==0 & e(sample)==1
replace meanasylum`i'_`z'_`x' = r(mean)
}
}
}
}

gen count = _n
#delimit ;
keep est* lb* ub* meanasylum* count;

#delimit cr
keep if count==1
#delimit ;
reshape long est1_f_ est2_f_ est3_f_ est4_f_ est5_f_ est6_f_ est7_f_ est8_f_ est9_f_ est10_f_ 
est1_m_ est2_m_ est3_m_ est4_m_ est5_m_ est6_m_ est7_m_ est8_m_ est9_m_ est10_m_
lb1_f_ lb2_f_ lb3_f_ lb4_f_ lb5_f_ lb6_f_ lb7_f_ lb8_f_ lb9_f_ lb10_f_ 
lb1_m_ lb2_m_ lb3_m_ lb4_m_ lb5_m_ lb6_m_ lb7_m_ lb8_m_ lb9_m_ lb10_m_
ub1_f_ ub2_f_ ub3_f_ ub4_f_ ub5_f_ ub6_f_ ub7_f_ ub8_f_ ub9_f_ ub10_f_ 
ub1_m_ ub2_m_ ub3_m_ ub4_m_ ub5_m_ ub6_m_ ub7_m_ ub8_m_ ub9_m_ ub10_m_
meanasylum1_f_ meanasylum2_f_ meanasylum3_f_ meanasylum4_f_ meanasylum5_f_ meanasylum6_f_ meanasylum7_f_ meanasylum8_f_ meanasylum9_f_ meanasylum10_f_ 
meanasylum1_m_ meanasylum2_m_ meanasylum3_m_ meanasylum4_m_ meanasylum5_m_ meanasylum6_m_ meanasylum7_m_ meanasylum8_m_ meanasylum9_m_ meanasylum10_m_, 
i(count) j(parents) string;

#delimit cr
drop count
gen count = _n
#delimit ;
reshape long est1_ est2_ est3_ est4_ est5_ est6_ est7_ est8_ est9_ est10_ 
lb1_ lb2_ lb3_ lb4_ lb5_ lb6_ lb7_ lb8_ lb9_ lb10_ 
ub1_ ub2_ ub3_ ub4_ ub5_ ub6_ ub7_ ub8_ ub9_ ub10_ 
meanasylum1_ meanasylum2_ meanasylum3_ meanasylum4_ meanasylum5_ meanasylum6_ meanasylum7_ meanasylum8_ meanasylum9_ meanasylum10_, 
i(count) j(gender) string;

#delimit cr
drop count
drop *6* *7* *8* *9* *10*
gen count = _n
#delimit ;
reshape long est lb ub meanasylum, 
i(count) j(quintile) string;

#delimit cr
gen quint = substr(quintile,1,1)
destring quint, replace

#delimit ;
twoway(scatter est quint if gender=="m_" & parents == "tp", mcolor(dkgreen))(rcap lb ub quint  if gender=="m_" & parents == "tp", color(dkgreen))
(line meanasylum quint  if gender=="m_" & parents == "tp", connect(l) yaxis(2) lpattern(dash) lcolor(orange_red)), yline(0, lpattern(solid) lcolor(gs8)) 
legend(off) xlabel(1 2 3 4 5) graphregion(fcolor(white)) xtitle("Parent Income Quintile") 
ytitle("Effect of Normal School County (Solid)", axis(1)) ytitle("Average Value in Asylum Counties (Dashed)", axis(2))
;
graph save working_men_tp.gph, replace;
graph export working_men_tp.pdf, replace;

#delimit ;
twoway(scatter est quint if gender=="m_" & parents == "sp", mcolor(dkgreen))(rcap lb ub quint  if gender=="m_" & parents == "sp", color(dkgreen))
(line meanasylum quint  if gender=="m_" & parents == "sp", connect(l) yaxis(2) lpattern(dash) lcolor(orange_red)), yline(0, lpattern(solid) lcolor(gs8)) 
legend(off) xlabel(1 2 3 4 5) graphregion(fcolor(white)) xtitle("Parent Income Quintile") 
ytitle("Effect of Normal School County (Solid)", axis(1)) ytitle("Average Value in Asylum Counties (Dashed)", axis(2))
;
graph save working_men_sp.gph, replace;
graph export working_men_sp.pdf, replace;

#delimit ;
twoway(scatter est quint if gender=="f_" & parents == "tp", mcolor(dkgreen))(rcap lb ub quint  if gender=="f_" & parents == "tp", color(dkgreen))
(line meanasylum quint  if gender=="f_" & parents == "tp", connect(l) yaxis(2) lpattern(dash) lcolor(orange_red)), yline(0, lpattern(solid) lcolor(gs8)) 
legend(off) xlabel(1 2 3 4 5) graphregion(fcolor(white)) xtitle("Parent Income Quintile") 
ytitle("Effect of Normal School County (Solid)", axis(1)) ytitle("Average Value in Asylum Counties (Dashed)", axis(2))
;
graph save working_women_tp.gph, replace;
graph export working_women_tp.pdf, replace;

#delimit ;
twoway(scatter est quint if gender=="f_" & parents == "sp", mcolor(dkgreen))(rcap lb ub quint  if gender=="f_" & parents == "sp", color(dkgreen))
(line meanasylum quint  if gender=="f_" & parents == "sp", connect(l) yaxis(2) lpattern(dash) lcolor(orange_red)), yline(0, lpattern(solid) lcolor(gs8)) 
legend(off) xlabel(1 2 3 4 5) graphregion(fcolor(white)) xtitle("Parent Income Quintile") 
ytitle("Effect of Normal School County (Solid)", axis(1)) ytitle("Average Value in Asylum Counties (Dashed)", axis(2))
;
graph save working_women_sp.gph, replace;
graph export working_women_sp.pdf, replace;
